{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation of the $\\mu$ opioid receptor with ligand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a complex build system as it has several components, the protein, a sodium ion, the ligand and of course the membrane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HTMD. 26-27 November 2015 HTMD workshop in Barcelona (fully booked)\n",
      "\n",
      "You are on the latest HTMD version (unpackaged).\n",
      "/tmp/tmp9sqlqd44\n"
     ]
    }
   ],
   "source": [
    "from htmd.console import *\n",
    "#get the files\n",
    "tmpdir = tempname()\n",
    "print(tmpdir)\n",
    "copytree(htmd.home()+'/data/mor/',tmpdir)\n",
    "chdir(tmpdir)\n",
    "path='./01_prepare/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['QM-min.pdb',\n",
       " '4dkl.pdb',\n",
       " 'sod.pdb',\n",
       " 'ff.rtf',\n",
       " 'ff.prm',\n",
       " 'membrane80by80C36.pdb']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glob('*')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found  segment between resid  65  and  263\n",
      "Found  segment between resid  270  and  352\n",
      "2015-11-16 10:32:43,649 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "Bond between A: [serial 3005 resid 140 resname CYS chain B segid P1]\n",
      "             B: [serial 3615 resid 217 resname CYS chain B segid P1]\n",
      "\n",
      "2015-11-16 10:32:44,027 - htmd.builder.builder - INFO - One disulfide bond was added\n",
      "2015-11-16 10:32:44,167 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2015-11-16 10:32:45,016 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "#Protein 4dkl is taken from opm\n",
    "\n",
    "topos  = ['top/top_all36_prot.rtf','top/top_all36_lipid.rtf', 'top/top_water_ions.rtf','ff.rtf']\n",
    "params = ['par/par_all36_prot.prm','par/par_all36_lipid.prm', 'par/par_water_ions.prm','ff.prm']\n",
    "prot = Molecule('4dkl.pdb')\n",
    "prot.filter('protein and noh and chain B or water within 5 of (chain B and protein)')\n",
    "pcenter = mean(prot.get('coords','protein'),axis=0)\n",
    "prot = segmentgaps(prot,'protein','P') \n",
    "unique(prot.get('segid'))\n",
    "prot = charmm.build(prot, topo=topos, param=params, outdir= path+'prot',ionize=False)\n",
    "# no need to change protonations\n",
    "#prot.view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2015-11-16 10:32:46,582 - htmd.builder.solvate - INFO - Using water pdb file at: /shared/sdoerr/Work/pyHTMD/htmd/builder/wat.pdb\n",
      "2015-11-16 10:32:46,933 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n",
      "Solvating: 100% (8/8) [############################################] eta 00:00 /\n",
      "2015-11-16 10:33:00,024 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "Bond between A: [serial 22800 resid 140 resname CYS chain B segid P1]\n",
      "             B: [serial 24036 resid 217 resname CYS chain B segid P1]\n",
      "\n",
      "2015-11-16 10:33:25,580 - htmd.builder.builder - INFO - One disulfide bond was added\n",
      "2015-11-16 10:33:25,699 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2015-11-16 10:33:26,929 - htmd.builder.charmm - INFO - Finished building.\n",
      "2015-11-16 10:33:27,854 - htmd.builder.ionize - INFO - Adding 14 anions + 0 cations for neutralizing and 70 ions for the given salt concentration.\n",
      "2015-11-16 10:33:28,122 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2015-11-16 10:33:28,122 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2015-11-16 10:33:28,123 - htmd.builder.ionize - INFO - Placing 84 ions.\n",
      "2015-11-16 10:33:51,898 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2015-11-16 10:34:13,701 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2015-11-16 10:34:15,003 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "#Add sodium in the receptor\n",
    "sod = Molecule('sod.pdb')\n",
    "sod.set('segid','S1')\n",
    "prot.append(sod)\n",
    "\n",
    "#Use a POPC membrane created with vmd and C36\n",
    "memb = Molecule('membrane80by80C36.pdb')\n",
    "mcenter = mean(memb.get('coords'),axis=0)\n",
    "memb.moveBy(pcenter-mcenter)\n",
    "mol = embed(prot,memb)\n",
    "\n",
    "#Add ligand, previously parametrized using gaamp\n",
    "lig = Molecule('QM-min.pdb') \n",
    "lig.set('segid','L')\n",
    "lcenter=mean(lig.get('coords'),axis=0)\n",
    "newlcenter=[random.uniform(-10, 10), random.uniform(-10, 10),  43 ]\n",
    "lig.rotateBy(uniformRandomRotation(), lcenter)\n",
    "lig.moveBy(newlcenter-lcenter)\n",
    "mol.append(lig) \n",
    "\n",
    "#Add water\n",
    "coo = mol.get('coords','lipids or protein')\n",
    "m = amin(coo,axis=0) + [0,0,-5]\n",
    "M = amax(coo,axis=0) + [0,0,20]\n",
    "mol = solvate(mol, minmax=vstack((m,M)))\n",
    "\n",
    "#Build\n",
    "mol = charmm.build(mol, topo=topos, param=params, outdir=path+'/build',saltconc=0.15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mol.view(sel='hetero',style='ball+stick', color='residueindex', hold=True)\n",
    "mol.view(viewer='vmd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.equilibration_v2 import Equilibration\n",
    "md = Equilibration()\n",
    "md.numsteps = 10000000\n",
    "md.temperature = 300\n",
    "md.fb_reference = 'protein and resid 293'\n",
    "md.fb_selection = 'segname L and noh'\n",
    "md.fb_box = [-25, 25, -25, 25, 43, 45]\n",
    "md.fb_k = 5\n",
    "md.acemd.celldimension = ' '.join(map(str,M-m))\n",
    "md.useconstantratio = 'on'\n",
    "md.write(path+'/build',path+'/equil')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mdx = LocalGPUQueue()\n",
    "mdx.submit(path+'/equil')\n",
    "mdx.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.production_v6 import Production\n",
    "md = Production()\n",
    "md.acemd.bincoordinates = 'output.coor'\n",
    "md.acemd.extendedsystem  = 'output.xsc'\n",
    "md.acemd.binvelocities=None\n",
    "md.acemd.binindex=None\n",
    "md.runtime = 50\n",
    "md.timeunits = 'ns'\n",
    "md.temperature = 300\n",
    "md.fb_reference = 'protein and resid 293'\n",
    "md.fb_selection = 'segname L and noh'\n",
    "md.fb_k = 5\n",
    "md.fb_box = [-25, 25, -25, 25, -10, 45]\n",
    "md.useconstantratio = 'on'\n",
    "md.write(path +'/equil','gen/s1')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}